Unsupervised Techniques

Principal Component Analysis (PCA)

Principal Component Analysis (PCA) is a popular dimensionality reduction technique used in Machine Learning applications. PCA condenses information from a large set of variables into fewer variables by applying some sort of transformation onto them. The transformation is applied in such a way that linearly correlated variables get transformed into uncorrelated variables. Correlation tells us that there is a redundancy of information and if this redundancy can be reduced, then information can be compressed. For example, if there are two variables in the variable set which are highly correlated, then, we are not gaining any extra information by retaining both the variables because one can be nearly expressed as the linear combination of the other. In such cases, PCA transfers the variance of the second variable onto the first variable by translation and rotation of original axes and projecting data onto new axes. The direction of projection is determined using eigenvalues and eigenvectors. So, the first few transformed features (termed as Principal Components) are rich in information, whereas the last features contain mostly noise with negligible information in them. This transferability allows us to retain the first few principal components, thus reducing the number of variables significantly with minimal loss of information.

Consider that you have a set of 2D points as it is shown in the figure above. Each dimension corresponds to a feature you are interested in. Here some could argue that the points are set in a random order. However, if you have a better look you will see that there is a linear pattern (indicated by the blue line) which is hard to dismiss. A key point of PCA is the Dimensionality Reduction. Dimensionality Reduction is the process of reducing the number of the dimensions of the given dataset. For example, in the above case it is possible to approximate the set of points to a single line and therefore, reduce the dimensionality of the given points from 2D to 1D.

Moreover, you could also see that the points vary the most along the blue line, more than they vary along the Feature 1 or Feature 2 axes. This means that if you know the position of a point along the blue line you have more information about the point than if you only knew where it was on Feature 1 axis or Feature 2 axis.

Hence, PCA allows us to find the direction along which our data varies the most. In fact, the result of running PCA on the set of points in the diagram consist of 2 vectors called eigenvectors which are the principal components of the data set.

The size of each eigenvector is encoded in the corresponding eigenvalue and indicates how much the data vary along the principal component. The beginning of the eigenvectors is the center of all points in the data set. Applying PCA to N-dimensional data set yields N N-dimensional eigenvectors, N eigenvalues and 1 N-dimensional center point.

PCA on Remote Sensing Imagery

Adjacent bands in a multi- or hyperspectral remotely sensed image are generally correlated. Multiband visible/NIR images of vegetated areas exhibit negative correlations between the NIR and visible red bands and positive correlations among the visible bands because the spectral reflectance characteristics of vegetation are such that as the vigour or greenness of the vegetation increases the red reflectance diminishes and the NIR reflectance increases. The presence of correlations among the bands of a multispectral image implies that there is redundancy in the data. Some information is being repeated. It is the repetition of information between the bands that is reflected in their intercorrelations.

Let's start by installing the necessary packages:

!pip install rasterio

Collecting rasterio
  Downloading rasterio-1.3.11-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (14 kB)
Collecting affine (from rasterio)
  Downloading affine-2.4.0-py3-none-any.whl.metadata (4.0 kB)
Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (24.2.0)
Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2024.8.30)
Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.7)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio) (0.7.2)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.26.4)
Collecting snuggs>=1.4.1 (from rasterio)
  Downloading snuggs-1.4.7-py3-none-any.whl.metadata (3.4 kB)
Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.1.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (71.0.4)
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.1.4)
Downloading rasterio-1.3.11-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (21.7 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 21.7/21.7 MB 24.6 MB/s eta 0:00:00
Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)
Downloading affine-2.4.0-py3-none-any.whl (15 kB)
Installing collected packages: snuggs, affine, rasterio
Successfully installed affine-2.4.0 rasterio-1.3.11 snuggs-1.4.7

import ee
ee.Authenticate()
ee.Initialize(project='my-project-1527255156007')

from google.colab import drive
drive.mount('/content/drive')

Mounted at /content/drive

import folium
from folium import plugins
from IPython.display import Image
import rasterio
import numpy as np
from matplotlib import pyplot as plt
import cv2
from matplotlib import cm
from collections import Counter
from skimage.transform import resize
from sklearn.cluster import KMeans

For this example we will use the image used in the last class:

path_barcelona = '/content/drive/MyDrive/Datasets_CV/Barcelona.tif'

with rasterio.open(path_barcelona) as src:
    im = src.read()

im = im.transpose([1,2,0])

img = im*2*255

img = img.astype('uint8')

red = img[:,:,2]
green = img[:,:,1]
blue = img[:,:,0]

And plot the RGB image:

rgb = np.dstack((red, green, blue))
plt.figure(figsize=[12,12])
plt.imshow(rgb)

<matplotlib.image.AxesImage at 0x7dbfacbdf1c0>

Before applying the PCA, it is necessary to resize the image by flattening the column and row dimensions into one:

pca_img = img.reshape((img.shape[0]*img.shape[1],img.shape[2]))

We import the SKlearn PCA:

from sklearn.decomposition import PCA
import pandas as pd

As we have 5 bands, let's set the number of components to 5:

pca = PCA(n_components = 5)

and then we perform the transformation:

principalComponents = pca.fit_transform(pca_img)

Let's see how much variance each component retains:

exp_variance = pca.explained_variance_ratio_

print(np.round(exp_variance,3))

[0.84  0.129 0.023 0.007 0.001]

Let's build a graph of components by variance and components by accumulated variance:

comp = np.arange(1,6)

df_variance  = pd.Series(exp_variance,name='explained variance')

df_cum_variance  = pd.Series(np.cumsum(exp_variance),name='cumulative explained variance')

df_comp = pd.Series(comp, name='Number of components')

df = pd.concat([df_variance, df_cum_variance, df_comp], axis=1)

import seaborn as sns
sns.set_theme(style="darkgrid")
from matplotlib import rcParams


rcParams['figure.figsize'] = 8,6

ax = sns.barplot(x="Number of components", y="explained variance", data=df)

ax = sns.lineplot(data=df, x="Number of components", y="cumulative explained variance")

Now we can plot each component separately:

principalComponents = principalComponents.reshape(img.shape[0],img.shape[1],img.shape[2])

fig = plt.figure(figsize = (30, 20))

for i in range(img.shape[2]):
    fig.add_subplot(1,img.shape[2], i+1)
    plt.imshow(principalComponents[:,:,i], cmap='inferno')
    plt.axis('off')
    plt.title(f'Band - '+ str(i+1))

We can also plot an image of the first three components as an RGB image:

PCA_nomr = cv2.normalize(principalComponents, np.zeros(principalComponents.shape),0,255 ,cv2.NORM_MINMAX)

PCA_nomr = PCA_nomr.astype('uint8')

plt.figure(figsize = (10, 10))
plt.imshow(PCA_nomr[:,:,0:3])
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

Clustering

Cluster analysis or clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense) to each other than to those in other groups (clusters).

Types of Clustering:

Broadly speaking, clustering can be divided into two subgroups :

• Hard Clustering: In hard clustering, each data point either belongs to a cluster completely or not. For example, in the above example each customer is put into one group out of the 10 groups.
• Soft Clustering: In soft clustering, instead of putting each data point into a separate cluster, a probability or likelihood of that data point to be in those clusters is assigned. For example, from the above scenario each costumer is assigned a probability to be in either of 10 clusters of the retail store

Types of clustering algorithms:

Since the task of clustering is subjective, the means that can be used for achieving this goal are plenty. Every methodology follows a different set of rules for defining the ‘similarity’ among data points. In fact, there are more than 100 clustering algorithms known. But few of the algorithms are used popularly, let’s look at them in detail:

• Connectivity models: As the name suggests, these models are based on the notion that the data points closer in data space exhibit more similarity to each other than the data points lying farther away. These models can follow two approaches. In the first approach, they start with classifying all data points into separate clusters & then aggregating them as the distance decreases. In the second approach, all data points are classified as a single cluster and then partitioned as the distance increases. Also, the choice of distance function is subjective. These models are very easy to interpret but lacks scalability for handling big datasets. Examples of these models are hierarchical clustering algorithm and its variants.
• Centroid models: These are iterative clustering algorithms in which the notion of similarity is derived by the closeness of a data point to the centroid of the clusters. K-Means clustering algorithm is a popular algorithm that falls into this category. In these models, the no. of clusters required at the end have to be mentioned beforehand, which makes it important to have prior knowledge of the dataset. These models run iteratively to find the local optima.
• Distribution models: These clustering models are based on the notion of how probable is it that all data points in the cluster belong to the same distribution (For example: Normal, Gaussian). These models often suffer from overfitting. A popular example of these models is Expectation-maximization algorithm which uses multivariate normal distributions.
• Density Models: These models search the data space for areas of varied density of data points in the data space. It isolates various different density regions and assign the data points within these regions in the same cluster. Popular examples of density models are DBSCAN and OPTICS.

Pixel clustering

Clustering can be used for unsupervised classification, which means that no training input is required, producing classes (i.e. clusters) that have no definition and consequently the user must assign a land cover label to each class.

KMeans on OpenCV:

OpenCV provides cv2.kmeans(samples, nclusters(K), criteria, attempts, flags) function for color clustering.

samples: It should be of np.float32 data type, and each feature should be put in a single column.

nclusters(K): Number of clusters required at the end

criteria: It is the iteration termination criteria. When this criterion is satisfied, the algorithm iteration stops. Actually, it should be a tuple of 3 parameters. They are ( type, max_iter, epsilon ):

Type of termination criteria. It has 3 flags as below:

cv.TERM_CRITERIA_EPS — stop the algorithm iteration if specified accuracy, epsilon, is reached. cv.TERM_CRITERIA_MAX_ITER — stop the algorithm after the specified number of iterations, max_iter. cv.TERM_CRITERIA_EPS + cv.TERM_CRITERIA_MAX_ITER — stop the iteration when any of the above condition is met.

attempts: Flag to specify the number of times the algorithm is executed using different initial labelings. The algorithm returns the labels that yield the best compactness. This compactness is returned as output.

flags: This flag is used to specify how initial centers are taken. Normally two flags are used for this: cv.KMEANS_PP_CENTERS and cv.KMEANS_RANDOM_CENTERS.

For this example we will use the rgb image from the previous example. We will also reduce the number of image dimensions from three to two:

vectorized = rgb.reshape((-1,3))

The next step is to convert to float32.

vectorized = np.float32(vectorized)

Let's apply the Kmeans from the openCV library. Let's split the image into 3 clusters:

K = 3
attempts=10
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)

ret,label,center=cv2.kmeans(vectorized,K,None,criteria,attempts,cv2.KMEANS_PP_CENTERS)

To reconstruct the clustered image, we take the resulting labels and resize them to the number of rows and columns of the original image:

label_img = label.reshape((rgb.shape[0:2]))

label_img = np.uint8(label_img)

Let's compare the result with the original image:

viridis = cm.get_cmap('viridis', K)

<ipython-input-38-689929f67468>:1: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  viridis = cm.get_cmap('viridis', K)

plt.figure(figsize=(14,14))
plt.subplot(1,2,1),
plt.imshow(rgb)
plt.title('Original Image')
plt.subplot(1,2,2),
plt.imshow(label_img,cmap=viridis)
plt.title('Clustering Image')
plt.show()

Color Quantization

Color Quantization is the process of reducing number of colors in an image. One reason to do so is to reduce the memory. Sometimes, some devices may have limitation such that it can produce only limited number of colors. In those cases also, color quantization is performed. Here we use k-means clustering for color quantization.

There is nothing new to be explained here. There are 3 features, say, R,G,B. So we need to reshape the image to an array of Mx3 size (M is number of pixels in image). And after the clustering, we apply centroid values (it is also R,G,B) to all pixels, such that resulting image will have specified number of colors. And again we need to reshape it back to the shape of original image.

center = np.uint8(center)

res = center[label.flatten()]
result_image = res.reshape((rgb.shape))

plt.figure(figsize=(14,14))
plt.subplot(1,2,1),
plt.imshow(rgb)
plt.title('Original Image')
plt.subplot(1,2,2),
plt.imshow(result_image)
plt.title('Result Image')
plt.show()

Clustering with Sklearn

to use Kmeans from SKlearn, let's use another image as an example.

path_vine = '/content/drive/MyDrive/Datasets_CV/DJI_0366.JPG'

with rasterio.open(path_vine) as src:
    img = src.read()

/usr/local/lib/python3.10/dist-packages/rasterio/__init__.py:304: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)

img.shape

(3, 3000, 4000)

img = img.transpose([1,2,0])

Let's plot a drone image over a viticulture:

plt.figure(figsize=[12,12])
plt.imshow(img)
plt.axis('off')

(-0.5, 3999.5, 2999.5, -0.5)

We followed practically the same procedure, reshape the image, setting the number of clusters to 3, obtaining the result of the clusters and reshape to the original value of the image:

img_vectorized = img.reshape((-1,3))

KMeans_pred = KMeans(n_clusters=3).fit_predict(img_vectorized)

/usr/local/lib/python3.10/dist-packages/sklearn/cluster/_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)

KMeans_pred = KMeans_pred.reshape((img.shape[0:2]))

As a result we have the image and three colors:

dark2 = cm.get_cmap('Dark2', 3)

<ipython-input-51-78bf2de087af>:1: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  dark2 = cm.get_cmap('Dark2', 3)

plt.figure(figsize=[12,12])
plt.imshow(KMeans_pred,cmap=dark2)
plt.axis('off')

(-0.5, 3999.5, 2999.5, -0.5)

After observing the resulting image, we can binarize the image by separating the soil from the plants and shadows:

binary_img = np.where(KMeans_pred==0,0,1)

plt.figure(figsize=[12,12])
plt.imshow(binary_img,cmap='Greys')
plt.axis('off')

(-0.5, 3999.5, 2999.5, -0.5)

Change Detection with PCA and Kmeans

Automatic change detection in images of a region acquired at different times is one the most interesting topics of image processing. Such images are known as multi temporal images. Change detection involves the analysis of two multi temporal satellite images to find any changes that might have occurred between the two time stamps. It is one of the major utilization of remote sensing and finds application in a wide range of tasks like defence inspections, deforestation assessment, land use analysis, disaster assessment and monitoring many other environmental/man-made changes.

We will be outlining an unsupervised method for change detection in this blog post. It involves the automatic analysis of the change data, i.e. the difference image, constructed using the multi temporal images. A difference image is the pixel-by-pixel subtraction of the 2 images. Eigen vectors of pixel blocks from the difference image will then be extracted by Principal Component Analysis (PCA). Subsequently, a feature vector is constructed for each pixel in the difference image by projecting that pixel’s neighbourhood onto the Eigen vectors. The feature vector space, which is the collection of the feature vectors for all the pixels, upon clustering by K-means algorithm gives us two clusters – one representing pixels belonging to the changed class, and other representing pixels belonging to the unchanged class. Each pixel will belong to either of the clusters and hence a change map can be generated. So, the steps towards implementing this application are:

1 - difference image generation and Eigen vector space (EVS)

2 - building the feature vector space (FVS)

3 - clustering of the feature vector space and change map

For this example, let's choose an area that has undergone some visible changes in recent years: Egypt's new capital.

AOI = ee.Geometry.Polygon(
        [[[31.442781230966744, 30.029336441993525],
          [31.442781230966744, 29.97641391019041],
          [31.538911602060494, 29.97641391019041],
          [31.538911602060494, 30.029336441993525]]])

We will select Sentinel 2 images from January 2016 and 2021

startDateviz_before = ee.Date.fromYMD(2016,1,1);
endDateviz_before = ee.Date.fromYMD(2016,1,31);
startDateviz_after = ee.Date.fromYMD(2021,1,1);
endDateviz_after = ee.Date.fromYMD(2021,1,31);

collectionviz_before = ee.ImageCollection("COPERNICUS/S2_HARMONIZED").filterDate(startDateviz_before,endDateviz_before).filterBounds(AOI).filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 1);
collectionviz_after = ee.ImageCollection("COPERNICUS/S2_HARMONIZED").filterDate(startDateviz_after,endDateviz_after).filterBounds(AOI).filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 1);

Then we select the median mosaic of the existing images and set the parameters for visualization:

S2_before = collectionviz_before.median().clip(AOI).divide(10000)
S2_after = collectionviz_after.median().clip(AOI).divide(10000)

vis_params = {'min': 0, 'max': 0.4, 'bands': ['B4', 'B3','B2']}

Here we define some basemaps to use in Folium:

basemaps = {
    'Google Satellite': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Satellite',
        overlay = True,
        control = True
    ),
    'Google Terrain': folium.TileLayer(
        tiles = 'https://mt1.google.com/vt/lyrs=p&x={x}&y={y}&z={z}',
        attr = 'Google',
        name = 'Google Terrain',
        overlay = True,
        control = True
    )}

The add_ee_layer function will make it easy to add vector and image layers to Folium:

def add_ee_layer(self, ee_object, vis_params, name):

    try:
        # display ee.Image()
        if isinstance(ee_object, ee.image.Image):
            map_id_dict = ee.Image(ee_object).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
            ).add_to(self)
        # display ee.ImageCollection()
        elif isinstance(ee_object, ee.imagecollection.ImageCollection):
            ee_object_new = ee_object.mosaic()
            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
            ).add_to(self)
        # display ee.Geometry()
        elif isinstance(ee_object, ee.geometry.Geometry):
            folium.GeoJson(
            data = ee_object.getInfo(),
            name = name,
            overlay = True,
            control = True
        ).add_to(self)
        # display ee.FeatureCollection()
        elif isinstance(ee_object, ee.featurecollection.FeatureCollection):
            ee_object_new = ee.Image().paint(ee_object, 0, 2)
            map_id_dict = ee.Image(ee_object_new).getMapId(vis_params)
            folium.raster_layers.TileLayer(
            tiles = map_id_dict['tile_fetcher'].url_format,
            attr = 'Google Earth Engine',
            name = name,
            overlay = True,
            control = True
        ).add_to(self)

    except:
        print("Could not display {}".format(name))

# Add EE drawing method to folium.
folium.Map.add_ee_layer = add_ee_layer

Now we can plot the two images and compare them:

my_map = folium.Map(location=AOI.centroid().coordinates().reverse().getInfo(), zoom_start=14)

# Add custom basemaps
basemaps['Google Terrain'].add_to(my_map)

# Add the elevation model to the map object.
my_map.add_ee_layer(S2_before, vis_params, 'New Cairo 2016')
my_map.add_ee_layer(S2_after, vis_params, 'New Cairo 2021')
my_map.add_child(folium.LayerControl())

# Display the map.
display(my_map)

Make this Notebook Trusted to load map: File -> Trust Notebook

Let's now download them to Drive and import using the rasterio:

image_before = S2_before.select(['B2','B3','B4','B5','B8']).reproject('EPSG:4326', None, 20)
image_after = S2_after.select(['B2','B3','B4','B5','B8']).reproject('EPSG:4326', None, 20)

task = ee.batch.Export.image.toDrive(image=image_before,
                                         crs='EPSG:4326',
                                         scale=20,
                                         fileFormat='GeoTIFF',
                                         description='New_Cairo_2016' ,
                                         maxPixels=1e13,
                                         folder='Datasets',
                                         region= AOI)

task.start()

task2 = ee.batch.Export.image.toDrive(image=image_after,
                                         crs='EPSG:4326',
                                         scale=20,
                                         fileFormat='GeoTIFF',
                                         description='New_Cairo_2021' ,
                                         maxPixels=1e13,
                                         folder='Datasets',
                                         region= AOI)

task2.start()

path_before = '/content/drive/MyDrive/Datasets_CV/New_Cairo_2016.tif'
path_after= '/content/drive/MyDrive/Datasets_CV/New_Cairo_2021.tif'

with rasterio.open(path_before) as src:
    im_before = src.read()

with rasterio.open(path_after) as src:
    im_after = src.read()

We will perform the same procedure we always do when importing new images:

im_before = im_before.transpose([1,2,0])
im_after = im_after.transpose([1,2,0])

im_before = im_before*2*255
im_after = im_after*2*255

im_before = im_before.astype('uint8')
im_after = im_after.astype('uint8')

Here, we are going to use the BGR to HSV conversion and use the intensity band to identify changes:

red_b = im_before[:,:,2]
green_b = im_before[:,:,1]
blue_b = im_before[:,:,0]

red_a = im_after[:,:,2]
green_a = im_after[:,:,1]
blue_a = im_after[:,:,0]

bgr_b = np.dstack((blue_b, green_b, red_b))
bgr_a = np.dstack((blue_a, green_a, red_a))

nir_before = im_before[:,:,3]
nir_after = im_after[:,:,3]

img_hsv_before = cv2.cvtColor(bgr_b, cv2.COLOR_BGR2HSV)
img_hsv_after = cv2.cvtColor(bgr_a, cv2.COLOR_BGR2HSV)

img_value_before = img_hsv_before[:,:,2]
img_value_after = img_hsv_after[:,:,2]

Now we can compare the images representing the intensity:

plt.figure(figsize=(18,14))
plt.subplot(1,2,1),
plt.imshow(img_value_before, cmap='inferno')
plt.axis('off')
plt.title('2016')
plt.subplot(1,2,2),
plt.imshow(img_value_after, cmap='inferno')
plt.title('2021')
plt.axis('off')
plt.show()

Let's create some helper functions to accomplish the task. The complete code is at the following github: https://github.com/abhijeet3922/Change-Detection-in-Satellite-Imagery/blob/master/PCAKmeans.py

def find_vector_set(diff_image, new_size):

    i = 0
    j = 0
    vector_set = np.zeros((int(new_size[0] * new_size[1] / 25), 25))

    while i < vector_set.shape[0]:
        while j < new_size[0]:
            k = 0
            while k < new_size[1]:
                block   = diff_image[j:j+5, k:k+5]
                #print(i,j,k,block.shape)
                feature = block.ravel()
                vector_set[i, :] = feature
                k = k + 5
            j = j + 5
        i = i + 1


    mean_vec   = np.mean(vector_set, axis = 0)
    vector_set = vector_set - mean_vec

    return vector_set, mean_vec

def find_FVS(EVS, diff_image, mean_vec, new):

    i = 2
    feature_vector_set = []

    while i < new[0] - 2:
        j = 2
        while j < new[1] - 2:
            block = diff_image[i-2:i+3, j-2:j+3]
            feature = block.flatten()
            feature_vector_set.append(feature)
            j = j+1
        i = i+1

    FVS = np.dot(feature_vector_set, EVS)
    FVS = FVS - mean_vec
    return FVS

def clustering(FVS, components, new):

    kmeans = KMeans(components, verbose = 0)
    kmeans.fit(FVS)
    output = kmeans.predict(FVS)
    count  = Counter(output)

    least_index = min(count, key = count.get)
    change_map  = np.reshape(output,(new[0] - 4, new[1] - 4))

    return least_index, change_map

First, let's resize the image to a size multiple of 5:

new_size = np.asarray(img_value_before.shape) / 5
new_size = new_size.astype(int) * 5

im_after = resize(img_value_after, (new_size))
im_before = resize(img_value_before, (new_size))

The conversion changed the range of values to 0 to 1. Let's return to the range from 0 to 255.

im_after = im_after*255
im_before = im_before*255

So we can create the difference image and then convert it to int16.

diff_image = abs(im_after - im_before)

diff_image = diff_image.astype('int16')

In this method, we take non-overlapping blocks of size 5 x 5 from the difference image and flatten them into row vectors:

vector_set, mean_vec = find_vector_set(diff_image, new_size)

PCA is then applied on this vector set to get the Eigen vector space. The Eigen vector space will be a 25 x 25 matrix; its each column is an Eigen vector of 25 dimensions.In Python, from sklearn.decomposition, we can simply import the PCA module and use it to perform PCA on vector_set variable to get the variable EVS:

pca = PCA()
pca.fit(vector_set)

PCA()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

EVS = pca.components_

Building the FVS involves again taking 5 x 5 blocks from the difference image, flattening them, and lastly projecting them onto the EVS, only this time, the blocks will be overlapping. A vector space (VS) is first made by constructing one vector for each pixel of the difference image such a way that one 5 x 5 block is actually a pixel’s 5 x 5 neighborhood. It is to be noted here that by this logic, 4 boundary rows and 4 boundary columns pixels won’t get any feature vectors since they won’t have a 5 x 5 neighborhood. (We can manage with this exclusion of these pixels, since it is safe to assume here that any changes occurring would be concentrated in the middle regions of the images, rather than the edges). So, we will have (m x n)- 8 feature vectors in the FVS, all 25 dimensional. Projecting the FVS to the 25 dimensional EVS simply means to perform the following matrix multiplication

FVS = find_FVS(EVS, diff_image, mean_vec, new_size)

The feature vectors for the pixels carry information whether the pixels have characteristics of a changed pixel or an unchanged one. Having constructed the feature vector space, we now need to cluster it so that the pixels can be grouped into two disjoint classes. We will be using the K-means algorithm to do that. Thus each pixel will get assigned to a cluster in such a way that the distance between the cluster’s mean vector and the pixel’s feature vector is the least. Each pixel gets a label from 1 to K, which denotes the cluster number that they belong to.

components = 3
least_index, change_map = clustering(FVS, components, new_size)

/usr/local/lib/python3.10/dist-packages/sklearn/cluster/_kmeans.py:1416: FutureWarning: The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning
  super()._check_params_vs_input(X, default_n_init=10)

Thus the argument components in clustering() will be 3. Remember, even though we have to do divide the pixels into 2 categories, we have chosen K = 3, instead of 2. Now how do we decide which of these clusters contains the pixels that belong to the changed class? It can be postulated that the cluster which contains the lowest number of pixels (denoted by variable least_index) is the cluster denoting the changed class, since the background remains more or less the same in satellite images and the changes occurred are comparatively less. Also, the mean of this cluster will be the highest. The reason behind the highest value of mean for that cluster is that the values of the difference image pixels in a region where some changes have occurred are higher than the values of pixels in the regions where there is no change.

Thus, in conclusion, the cluster with the lowest number of pixels, and also the highest mean is the cluster belonging to the changed class.

With this information, we will now build a change map – a binary image to show the output of change detection. We have chosen to keep the background black and will show the changes in white, i.e., intensity value of those pixels will be 255. You can do the reverse as well.

change_map[change_map == least_index] = 255
change_map[change_map != 255] = 0
change_map = change_map.astype(np.uint8)

Finally, let's erode the image to filter out some noise:

kernel = np.asarray(((0,0,1,0,0),
                     (0,1,1,1,0),
                     (1,1,1,1,1),
                     (0,1,1,1,0),
                     (0,0,1,0,0)), dtype=np.uint8)

cleanChangeMap = cv2.erode(change_map,kernel)

Finally, we can present the difference image and the change map of the study region:

plt.figure(figsize=(18,14))
plt.subplot(1,2,1),
plt.imshow(diff_image, cmap='Greys_r')
plt.axis('off')
plt.title('Difference Image')
plt.subplot(1,2,2),
plt.imshow(cleanChangeMap, cmap='Greys_r')
plt.title('Change Map')
plt.axis('off')
plt.show()

Unsupervised Segmentation

In digital image processing and computer vision, image segmentation is the process of partitioning a digital image into multiple segments (sets of pixels, also known as image objects). The goal of segmentation is to simplify and/or change the representation of an image into something that is more meaningful and easier to analyze. Image segmentation is typically used to locate objects and boundaries (lines, curves, etc.) in images. More precisely, image segmentation is the process of assigning a label to every pixel in an image such that pixels with the same label share certain characteristics.

The result of image segmentation is a set of segments that collectively cover the entire image, or a set of contours extracted from the image (see edge detection). Each of the pixels in a region are similar with respect to some characteristic or computed property, such as color, intensity, or texture. Adjacent regions are significantly different color respect to the same characteristic(s). When applied to a stack of images, typical in medical imaging, the resulting contours after image segmentation can be used to create 3D reconstructions with the help of interpolation algorithms like marching cubes.

Unsupervised segmentation requires no prior knowledge. Consider an image that is so large that it is not feasible to consider all pixels simultaneously. So in such cases, Unsupervised segmentation can breakdown the image into several sub-regions, so instead of millions of pixels, you have tens to hundreds of regions.

from skimage.segmentation import felzenszwalb, slic, quickshift, watershed
from skimage.segmentation import mark_boundaries
from skimage.util import img_as_float
from skimage.filters import sobel
from skimage.color import rgb2gray
from skimage import color

Let's prepare our RGB image, converting to float:

img = img_as_float(rgb)

Felzenszwalb’s efficient graph based segmentation

This fast 2D image segmentation algorithm, proposed in 1 is popular in the computer vision community. The algorithm has a single scale parameter that influences the segment size. The actual size and number of segments can vary greatly, depending on local contrast.

segments_fz = felzenszwalb(img, scale=100, sigma=0.5, min_size=50)

Let's see how the segments were divided:

vals = np.linspace(0,1,250)
np.random.shuffle(vals)
cmap = plt.cm.colors.ListedColormap(plt.cm.jet(vals))
plt.figure(figsize=[12,12])
plt.imshow(segments_fz, cmap=cmap )
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

We can also overlap the outlines over the original image:

plt.figure(figsize=[12,12])
plt.imshow(mark_boundaries(img, segments_fz))
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

SLIC - K-Means based image segmentation

This algorithm simply performs K-means in the 5d space of color information and image location and is therefore closely related to quickshift. As the clustering method is simpler, it is very efficient. It is essential for this algorithm to work in Lab color space to obtain good results. The algorithm quickly gained momentum and is now widely used. See 3 for details. The compactness parameter trades off color-similarity and proximity, as in the case of Quickshift, while n_segments chooses the number of centers for kmeans.

segments_slic = slic(img, n_segments=250, compactness=10, sigma=1, start_label=1)

plt.figure(figsize=[12,12])
plt.imshow(mark_boundaries(img, segments_slic))
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

Let's create an image representing the average value of the RGB channels per segment:

LabelsAve = color.label2rgb(segments_slic, img, kind='avg')
plt.figure(figsize=[12,12])
plt.imshow(LabelsAve)
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

Quickshift image segmentation

Quickshift is a relatively recent 2D image segmentation algorithm, based on an approximation of kernelized mean-shift. Therefore it belongs to the family of local mode-seeking algorithms and is applied to the 5D space consisting of color information and image location 2.

One of the benefits of quickshift is that it actually computes a hierarchical segmentation on multiple scales simultaneously.

Quickshift has two main parameters: sigma controls the scale of the local density approximation, max_dist selects a level in the hierarchical segmentation that is produced. There is also a trade-off between distance in color-space and distance in image-space, given by ratio.

segments_quick = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)

plt.figure(figsize=[12,12])
plt.imshow(mark_boundaries(img, segments_quick))
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

Compact watershed segmentation of gradient images

Instead of taking a color image as input, watershed requires a grayscale gradient image, where bright pixels denote a boundary between regions. The algorithm views the image as a landscape, with bright pixels forming high peaks. This landscape is then flooded from the given markers, until separate flood basins meet at the peaks. Each distinct basin then forms a different image segment. 4

As with SLIC, there is an additional compactness argument that makes it harder for markers to flood faraway pixels. This makes the watershed regions more regularly shaped

gradient = sobel(rgb2gray(img))
segments_watershed = watershed(gradient, markers=250, compactness=0.001)

plt.figure(figsize=[12,12])
plt.imshow(mark_boundaries(img, segments_watershed))
plt.axis('off')

(-0.5, 816.5, 512.5, -0.5)

You can also convert the segments to vectors and save it as a SHAPEFILE file for example:

segments_watershed = segments_watershed.astype('uint16')

Let's import the "features" function from the rasterio and the geopandas and shapely libraries

from rasterio import features
from shapely.geometry import shape
import geopandas as gpd
import shapely

First we perform the image to vector conversion:

src_barc = rasterio.open(path_barcelona)

find_shapes = features.shapes(segments_watershed, transform=src_barc.transform)

Then we get the coordinates of each vector:

polygons = [shapely.geometry.Polygon(shape[0]["coordinates"][0]) for shape in find_shapes]

We created a geodataframe and set the crs:

crs = {'init': 'epsg:4326'}
my_gdf = gpd.GeoDataFrame(crs=crs, geometry=polygons)

/usr/local/lib/python3.10/dist-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  in_crs_string = _prepare_from_proj_string(in_crs_string)

And finally we save the file on the drive:

my_gdf.to_file("/content/Segments.shp", driver='ESRI Shapefile')

Apply maskSLIC

The maskSLIC method is an extension of the SLIC method for the generation of superpixels in a region of interest. maskSLIC is able to overcome border problems that affects SLIC method, particularely in case of irregular mask.

Let's use the image, the image of a viticulture and the mask of the region with valid pixels:

path_vine2 = '/content/drive/MyDrive/Datasets_CV/Drone_img.tif'

src = rasterio.open(path_vine2)
img2 = src.read()

img2 = img2.transpose([1,2,0])

img2 = img2.astype('uint8')

<ipython-input-107-3b19470436d4>:1: RuntimeWarning: invalid value encountered in cast
  img2 = img2.astype('uint8')

plt.figure(figsize=[12,12])
plt.imshow(img2)
plt.axis('off')

(-0.5, 3793.5, 2270.5, -0.5)

plt.figure(figsize=[12,12])
plt.imshow(img2[:,:,3])
plt.axis('off')

(-0.5, 3793.5, 2270.5, -0.5)

We apply the SLIC setting the mask and segmentation occurs only in the valid region:

mask = img2[:,:,3]

img2 = img2[:,:,0:3]

img2_float = img_as_float(img2)

m_slic = slic(img2_float, n_segments=500, mask=mask, start_label=1)

plt.figure(figsize=[12,12])
plt.imshow(mark_boundaries(img2, m_slic))
plt.contour(mask, colors='red', linewidths=1)
plt.axis('off')

(-0.5, 3793.5, 2270.5, -0.5)

References:

https://www.analyticsvidhya.com/blog/2016/11/an-introduction-to-clustering-and-different-methods-of-clustering/

https://www.kdnuggets.com/2019/08/introduction-image-segmentation-k-means-clustering.html

https://docs.opencv.org/4.5.2/d1/d5c/tutorial_py_kmeans_opencv.html

https://en.wikipedia.org/wiki/Image_segmentation

https://docs.opencv.org/master/d3/db4/tutorial_py_watershed.html

https://towardsdatascience.com/image-segmentation-using-pythons-scikit-image-module-533a61ecc980

https://scikit-image.org/docs/stable/auto_examples/segmentation/plot_segmentations.html

https://scikit-image.org/docs/dev/auto_examples/segmentation/plot_mask_slic.html

https://towardsdatascience.com/principal-component-analysis-in-depth-understanding-through-image-visualization-892922f77d9f

https://docs.opencv.org/3.4/d1/dee/tutorial_introduction_to_pca.html

https://towardsdatascience.com/comprehensive-guide-to-satellite-imagery-analysis-using-python-1b4153bad2a

https://github.com/abhijeet3922/Change-Detection-in-Satellite-Imagery/blob/master/PCAKmeans.py

https://appliedmachinelearning.blog/2017/11/25/unsupervised-changed-detection-in-multi-temporal-satellite-images-using-pca-k-means-python-code/